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EARTHQUAKE SIZE DISTRIBUTION: POWER-LAW WITH EXPONENT = \ ? 
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Abstract. We propose that the widely observed and universal Gutenberg- Richter relation 
is a mathematical consequence of the critical branching nature of earthquake process in a 
brittle fracture environment. These arguments, though preliminary, are confirmed by recent 
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systematic and random errors in determining earthquake size, especially its seismic moment. 
These effects increase the estimate of the parameter /3 of the power-law distribution of 
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earthquake sizes. In particular, we find that estimated /3-values may be inflated by 1-3% 

(N 

because relative moment uncertainties decrease with increasing earthquake size. Moreover, 
earthquake clustering greatly influences the /3-parameter. If clusters (aftershock sequences) 
are taken as the entity to be studied, then the exponent value for their size distribution would 



decrease by 5-10%. The complexity of any earthquake source also inflates the estimated /3- 
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value by at least 3-7%. The centroid depth distribution also should influence the /3-value, 
an approximate calculation suggests that the exponent value may be increased by 2-6%. 
Taking all these effects into account, we propose that the recently obtained /3-value of 0.63 
could be reduced to about 0.52-0.56: near the universal constant value (1/2) predicted 
by theoretical arguments. We also consider possible consequences of the universal /3-value 
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and its relevance for theoretical and practical understanding of earthquake occurrence in 
various tectonic and Earth structure environments. Using comparative crystal deformation 
results may help us understand the generation of seismic tremors and slow earthquakes and 
illuminate the transition from brittle fracture to plastic flow. 

Short running title: EARTHQUAKE SIZE DISTRIBUTION 

Key words: Gutenberg-Richter relation; Corner moment; Tapered Pareto distribution; 
Scalar and tensor seismic moment; Universality of earthquake size distribution; Random 
walk; 3-D random rotation; Earthquake depth distribution; Seismic tremors; Transition 
from brittle to plastic deformation. 

1 Introduction 

Earthquake size distribution is usually described by the Gutenberg-Richter (G-R) magnitude- 
frequency relation (see Section I2~2"j) . The G-R distribution can be transformed into the Pareto 
(power-law) distribution for a scalar seismic moment M with the exponent = 6/1.5, where 
b is the parameter of the G-R law. 

Theoretical analysis of earthquake occurrence (Vere- Jones, 1976, 1977) suggests that, 
given its branching nature, the exponent /3 of earthquake size distribution should be identical 
to 1/2. Properties of the critical branching process explain this result: the total number of 
events (individuals) in such a process asymptotically is distributed according to a power-law 
with exponent 0.5 (Otter, 1949; Harris, 1963, Ch. 1.13). Such distributions, obtained by 
simulations, are shown as upper curves in Figs. [7] and [8] below. 

The same values of power-law exponents are derived for percolation and self-organized 
critcality (SOC) processes in a high- dimensional space (see discussion by Kagan, 1991a, 
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p. 132). Similar values of exponents are obtained by theoretical arguments and simulations 
for dislocation avalanches in crystalline materials (Zaiser, 2006, and references therein). 

However, almost all the (3 or b measurements in earthquake catalogs result in estimates 
larger than | (or 0.75 for the b- value). These estimates exhibit large variations in different 
regions, tectonic zones, etc. A search of the ISI database (Thomson Reuters Scientific) 
indicates that over the last five years, more than three papers on the 6-value or the G-R 
relation were published monthly. The similar rate of publications of the papers discussing 
the magnitude-frequency relation can be observed in the previous 20-30 years. In almost all 
articles the variation of the 6-value is attributed to different tectonics, rock stress, etc. 

The following reasons for variability in the measured b- and /3-values can be proposed: 

• 1. Inappropriate usage of magnitude scales other than moment magnitude: only the 
moment magnitude should be studied. Regular earthquake magnitudes have significant 
systematic and random errors (Kagan, 1999, 2003), making them inappropriate for rigorous 
statistical, quantitative investigation. 

• 2. The maximum or corner moment M c (see Section 12.21) needs careful consideration 
(Kagan, 1991a; 2002a; Bird and Kagan, 2004). If M c close to the magnitude /moment thresh- 
old - that is, the smallest magnitude/moment above which the catalog can be considered 
to be complete - and the earthquake size distribution is approximated by a plain G-R dis- 
tribution, the magnitude-frequency curve would shift downwards, and the b- or /3-estimates 
would be strongly biased upwards. In such a case to avoid the bias it is necessary to apply 
a two-parameter relation, which includes the maximum or corner magnitude. 

• 3. Mixing populations of earthquakes from tectonic settings that have different corner 
magnitudes: when earthquake populations with a varying corner moment M c are placed in 
the study samples, a false increase in the fi- values results. Mixing populations with different 
M c may even yield a seemingly linear curve in a log-log plot, in such two-parameter 
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approximation of the magnitude-frequency relation would fail to avoid bias. For example, due 
to a significant corner moment variability, Kagan (2002a, Table 5, Section 5.2.3) determined 
the (3- values on the order 0.8-1.1 for mid-ocean earthquakes. Similarly and apparently for 
the same reason, excessively large 6-values for oceanic earthquakes were obtained by Okal 
and Romanowicz (1994) and Schorlemmer et al. (2005). 

• 4. Relative seismic moment errors increase with decreases in earthquake size, resulting 
in a spurious (3 increase. 

• 5. The object of study should be earthquake sequences, not individual earthquakes: 
theoretical estimates, discussed above, are relevant for earthquake sequences, not individual 
events. Hence (3- values need to be corrected for this effect. 

• 6. An earthquake is a tensor; its size, as given in moment tensor solutions, is a tensor 
sum of earthquake subevents. If these subevents have different focal mechanisms, their 
tensor sum would be smaller than that for scalar moments of subevents. Thus, even if the 
number of elementary earthquake events were distributed according to the power-law with 
the exponent /3 = |, the distribution of earthquake size, as represented by the tensor sum, 
may have a larger exponent value. This follows from stochastic complexity in the source. If, 
for example, a source consists of positive and negative random n subsources, its size would 
be proportional to y/n: the size would be similar to that at the end of a Brownian random 
walk. 

• 7. The corner moment M c is likely to change significantly with the depth for shallow 
earthquakes (Kagan, 1999, Section 1; Kagan, 2002a). For shallow seismicity, earthquake 
populations with various depths are usually added up, thus j3 determinations are biased. 
Earthquake catalogs with only the hypocenter information cannot be used to investigate 
this effect because hypocenters are often located at a lower or upper depth boundary of a 
fault rupture zone. On the other hand, the degree of accuracy of centroid depths in moment 
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tensor catalogs is presently too low to study rigorously the M c depth variation. 

On the basis of statistical analysis of several earthquake catalogs and some theoretical 
observations, Kagan (1991a, pp. 129, 132-3) conjectured that the /3-value is a universal 
constant (1/2). Three additional arguments can now be added to strengthen this hypothesis: 
(a) more careful measurements of the ft parameter in modern earthquake catalogs (mainly 
the CMT catalog) suggest that the ft is universal; (b) recent advances in space geodesy and 
quantitative plate tectonics allowed for detailed calculation of tectonic deformation rate. By 
comparing tectonic and seismic moment rates, we can calculate the upper bound, M c , for 
earthquake moment distribution; and (c) investigation of dislocation avalanches in crystals 
indicated that their size distribution is a power-law whose exponent has a universal value. 

Several previous investigations (Kagan, 1999; Bird and Kagan, 2004; Boettcher and 
Jordan, 2004) suggested that the ft has a universal value on the order of 0.63-0.67. Moreover, 
statistical analysis of earthquake size distribution and comparison of seismic and tectonic 
deformation rates allowed us to evaluate the corner moment (M c ) value for several tectonic 
regions. We conclude that apparent change in the /3-values is due mainly to M c variability. 

However, the ft universality model can be challenged. One can argue that more careful 
measurements may reveal statistically significant variations in the ft exponent. The aim 
here is to show that the ft- value is a universal constant. We consider various systematic and 
random effects whose influence would confirm this conjecture. 

Recent experimental and theoretical investigations have demonstrated that crystal plas- 
ticity is characterized by large intrinsic spatio-temporal fluctuations with scale-invariant 
characteristics. In other words, deformation proceeds through intermittent bursts with 
power-law size distributions (Zaiser, 2006; Dahmen et al, 2009). In particular, Csikor et 
al. (2007), Dimiduk et al. (2006), Weiss and Marsan (2003) studied dislocation avalanches 
(micro- and nano-earthquakes) in ice and other crystals. They consistently obtain the power- 
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law size distribution with the probability density exponent close to 1.5. This would corre- 
spond to our exponent 1 + (3. The power-size distribution in a single ice crystal extends over 
six decades of magnitude (Miguel et al, 2001). 

Comparing these dislocation avalanche measurements with seismological observations 
leads to some problems. The most accurate measurements of earthquake size are for the 
seismic moment. Experimental laboratory observations provide the energy of acoustic emis- 
sion bursts or strain step measurements. Zaiser (2006, pp. 212, 223) argues that both these 
measurements are approximately equivalent and yield similar values for the distribution 
density power-law exponent (1.5-1.6). 

Earthquake energy has the same distribution as the seismic moment. Kanamori (1977) 
cites the following relation between energy released by earthquakes, E, and their magnitude, 
m: log 10 E = 1.5 m + 11.8. Because the moment and magnitudes have a similar dependence 
(Kagan, 2002a), the exponents for energy and moment power-law distributions should be 
identical. 

Another problem of comparison is that obtained statistical distributions of dislocation 
avalanches are not processed by appropriate statistical techniques. Thus, the obtained values 
of the exponents may be biased. For example, Zaiser et al. (2008) use the least-square fit 
to calculate the exponent. Richeton et al. (2005) apply the Levenberg-Marquardt algorithm 
for this purpose. These methods are appropriate for fitting regression curves, but they 
are not the most statistically efficient techniques for a parameter evaluation in statistical 
distribution. They may yield biased estimates of an exponent parameter and its uncertainty 
(Vere- Jones, 1988; Clauset et al., 2009). 

In Section [2]below we briefly review the earthquake size distribution formulas and consider 
the results of the (3- value evaluation. Section [3] discusses systematic and random effects in 
earthquake size determination (items 4-7 above). In Section H] we summarize the results 
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for determining /3 and compare our conclusions with studies of dislocation avalanches in 
brittle and plastic crystalline solids. We also discuss possible consequences for interpreting 
geophysical observations of regular and slow earthquakes and seismic tremors in brittle crust 
and the upper mantle. 

2 Catalog analysis and earthquake size distribution 

Only shallow earthquakes (depth 0-70 km) will be investigated in this work, because more 
data is available on them. Additionally, the seismic efficiency coefficient, or proportion 
of tectonic deformation released by such events, is close to 1.0 (Bird and Kagan, 2004). 
Geometry of deep earthquake faults is much less known than that for shallow seismicity, and 
may be more complex. For deeper earth layers (depth > 70 km) only a small part (less than 
5%) of tectonic motion is released by seismicity (Kagan, 1999; Frohlich and Wetzel, 2007). 
Therefore, earthquake rupture in a brittle crust would be better modeled by the critical 
branching process theory. 

As mentioned in the Introduction, regular magnitude measurements are subject to many 
random and systematic errors (Kagan, 2003). Kagan (1999, pp. 557-8) studied correlation 
between 6-value estimates based on m b and M s magnitudes and found that the correlation 
coefficient is low (0.1-0.2). These coefficient values seem to indicate that b- value variations 
are not caused by regional tectonic or physical factors. Therefore, we investigate earthquake 
size distribution here using only seismic moment tensor measurements. The most complete, 
extensive, and accurate catalog of tensor solutions is the CMT dataset (Kagan, 2003). 
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2.1 CMT earthquake catalog 

We studied earthquake distributions and clustering for the global CMT catalog of moment 
tensor inversions compiled by the CMT group (Ekstrom et ai, 2005; Ekstrom, 2007). The 
present catalog contains more than 30,000 earthquake entries for the period 1977/1/1 to 
2008/12/31. Earthquake size is characterized by a scalar seismic moment M. The moment 
magnitude can be calculated from the seismic moment (measured in Nm) value as 

m w = (2/3) • log 10 M - 6.0 . (1) 

The magnitude threshold for the 1977-2008 catalog is m t = 5.8; for the 1982-2008 catalog it 
is m t = 5.6 (Kagan, 2003). An earthquake catalog is considered reasonably complete at and 
above the magnitude threshold (or the corresponding seismic moment, Eq. [1]). Since we use 
only the moment magnitude in this work, the subscript is usually omitted. 

2.2 Seismic moment distribution 

The distribution of earthquake size is usually described by the classical G-R (Gutenberg and 
Richter, 1944) magnitude-frequency relation 

log 10 N(m) = a t — b(m — m t ) for m t < m, (2) 

where N(m) is the number of earthquakes with magnitude > m, and a t and b are parameters: 
a t is the logarithm of the number of earthquakes with m > m t (the seismic activity or 
earthquake productivity level) and b describes the relation between numbers of small and 
large earthquakes. The G-R laws seems to apply to earthquakes as small as m = —1.3 
(Boettcher et ai, 2009), with a rupture length on the order of 0.5 m, and as large as the 
2004 m9.2 Sumatra, with about a 1200 km rupture. 



The G-R relation (jSJ) can be transformed into the Pareto (power-law) distribution for the 
scalar seismic moment M. The distribution in a probability density form is (Kagan, 2002a) 

f(M) = fiMfM- 1 -? for M, < M : (3) 

in a survivor function (1 — cumulative distribution) form it is 

F(M) = (Mt/Mf for M t < M, (4) 

where (3 is the index parameter of the distribution, and b — | /3. 

The tapered G-R (TGR) relation has an exponential taper applied to the number of 
events with a large seismic moment. Its probability density function is 

U1 Mf Af-^expf^ ~ M ) for M t < M < oo . 



f(M) 



M, 



c 



(5) 



Here M c is the parameter controlling the distribution in the upper ranges of M ('the corner 
moment'). The survivor function is 

F(M) = {Mt/Mf exp ( Mt ~ M \ for M t < M < oo . (6) 



M c 

Equations (EJ) and <$5§ are equivalent to (J2D and PJ, respectively, if M c — > cxo. 

Fig. [1] displays cumulative distribution (survivor function) for the scalar seismic moment 
of shallow earthquakes in the CMT catalog for 1977-2008. The curves display a scale- 
invariant (Pareto) segment (linear in the log-log plot) for small and moderate values of the 
moment magnitude m. But for large m, the curve clearly bends downward (see Fig. 2 by 
Kagan et ai, 2010 for comparison). 

In a regular /3 estimation we construct a likelihood map in [(3 x m c ] space (Fig. and 
find its maximum (compare Figs. 6 and 7 by Bird and Kagan, 2004, or Fig. 3 by Kagan et 
ai, 2010). The 95% confidence level corresponds to the contour value 0.0. Only the lower 
confidence limit for the corner magnitude m c ~ 8.3 can be determined. The upper limit is 
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not denned in the map. The tapered Pareto distribution with the parameters corresponding 
to the maximum likelihood function, is shown in Fig. [T]by the dashed line. Sometimes these 
maps exhibit a significant correlation between the estimates of parameters, see Figs. 6 G,H 
and 7 B,E by Bird and Kagan (2004), complicating M c evaluation and application. Since our 
map does not display the interdependence between these two parameter estimates, simpler 
procedures for determining (3 can be used (Kagan, 2002a). 

In Fig. [3] we show the normalized difference between the observed magnitude-frequency 
relation in Fig. [1] and its approximation by the tapered Pareto or the tapered Gutenberg- 
Richter (TGR) relation. We assume that the earthquake numbers in a frequency plot follow 
the Poisson distribution. Kagan (2010) demonstrates that the temporal distribution of large 
earthquakes approaches the Poisson law. Kagan et ai, (2010, Table 2) show that in the 
CMT catalog due to its high moment threshold the aftershock numbers are less than 25% 
of the total. For large number of events its standard deviation (a) is a square root of the 
number. To normalize, we divide the difference by 1.96 x a; hence the ±1 ordinate value 
would correspond to the 95% confidence level. 

Fig. 4 demonstrates a possible source of bias in determining magnitude: earthquakes of 
different size have their size estimated differently The moment inversion in the CMT catalog 
is carried out using some combination of three types of waves: body, surface, and mantle 
waves (Dziewonski et ai, 1981; Dziewonski and Woodhouse, 1983; Ekstrom, 2007). The 
mantle low-frequency (period 135 s) waves are used mostly for larger earthquakes. Higher- 
frequency body and surface waves are used for more moderate events (see also Fig. 2 by 
Kagan, 2003). For shallow earthquakes, Fig. 0] shows the fraction of each wave used in the 
solution for each moment magnitude. Earthquake sizes are estimated by some combination 
of the different wave types, so the sum of three points in any given magnitude bin may exceed 
unity. For example, up to m = 7.5 almost 90% of earthquakes have body waves implemented 
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in the solution. After that the inversion is based mostly on mantle waves. Therefore, we 
may posit that the statistically significant deviation for the tapered Pareto distribution at 
m > 7.8 (see Fig. [3]) may be caused by transition from the body- to mantle-wave estimate. 
It is possible, of course, that the change in both plots (Figs. [3]and|l]) at about m = 7.5 is a 
coincidence. Presently we lack data to thoroughly test this hypothesis. 



3 Systematic and random effects in determining earth- 
quake size 

This section considers in detail four sources of magnitude/moment bias mentioned in the 
Introduction (items 4-7). 

3.1 Scalar seismic moment errors 

There is a bias in evaluating the G-R parameters due to random errors in determining mag- 
nitude (Molchan and Podgaetskaya, 1973, p. 47, Eq. 7; Tinti and Mulargia, 1985, p. 1690). 
This magnitude error is assumed to be symmetric and Gaussian, the error analysis by Kagan 
(2003, see Figs. 12-15) seems to confirm this. Generally, as long as the moment or amplitude 
errors are relatively small, one should expect them to be Gaussian as they are the result of 
summing up many independent random uncertainties. 

The magnitude error causes a shift of the a estimate in (T5]) toward larger values: given 
the approximate symmetry of the error distribution, more weak earthquakes have their 
magnitude increased than vice versa: 

b 2 a 2 log 10 

Westim w corr 2 ' 

where a m is a standard error in the magnitude estimate, a corr is the corrected (true) a- value, 
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Oestim is the measured a- value. Otherwise, the bias may manifest as a right-hand (horizontal) 
magnitude shift in the G-R curve 



where m corr is the corrected magnitude, m est i m is the measured magnitude. For b = 1 the 
shift is of the same amplitude in (J7J) and (jHJ). 

If the magnitude errors do not depend on the magnitude, this error does not practically 
influence the estimated b- value (Tinti and Mulargia, 1985). However, if function of 

m, the b estimates would be affected by magnitude errors. 

Rhoades (1996, 1997) derived the theoretical estimates of a bias in the b- value in such a 
case. Rhoades (1996) and Rhoades and Dowrick (2000) studied the influence of magnitude 
errors on b- value estimates and provided some approximate estimates for the 6-bias from such 
errors. In these evaluations they assumed that magnitude errors increase as the magnitude 
itself increases: for the j-th measurement of magnitude Oj = 0.1 ( 1 + Uj rrij ), where Uj is a 
random number uniformly distributed on the interval [0, 1], the magnitude threshold is 3.95 
and the real 6-value is 1.0. They found that the 'measured' 6-value decreases by about 4% 
compared to its true value. 

The CMT catalog supplies inversion errors for tensor moment components (Dziewonski 
et ai, 1981; Dziewonski and Woodhouse, 1983). It is important to distinguish magnitude 
and moment random errors. Magnitude is defined as a logarithm of appropriately scaled 
and averaged amplitude (A) of specific seismic waves: m = log 10 A + C, where C is a 
coefficient. Therefore, for small amplitude uncertainties the standard error in the amplitude 
measurements (a a) is related to the magnitude error (<r m ) as 



ba 2 m log 10 



(8) 



m, 



corr 



2 



dm 




(9) 



ni 



oc 
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i.e., the magnitude error is proportional to a relative amplitude error. A similar relation is 
valid for the scalar seismic moment (see Eqs. ITUl and [T3l below) . 

Kagan (2002a) measured the scalar seismic moment errors in the CMT catalog and found 
that the relative moment errors actually decrease with the earthquake size, implying that 
the magnitude errors should also decrease. We define the relative moment error, e, as the 
ratio of the error tensor norm to the moment tensor norm 



where Eij and My are standard error and moment tensor components, respectively. The 
distribution of e for the CMT catalog 1977-2008 is shown in Fig. [5] (compare Fig. 5 in Kagan, 
2002a). Since the influence of the magnitude threshold value is insignificant for determining 
the relative moment error, we use m t = 5.6 in this plot. We calculate two regression lines 
approximating the dependence of the errors on the magnitude: the linear and quadratic 
curves 



We use (m — 6) instead of m as an argument, so that the cq- value would have a clear intuitive 
meaning. 

In the diagram (Fig. [5]) the coefficient of correlation between e and the magnitude is 
—0.47, indicating that relative moment errors decrease with the increase of m. Residual 
regression errors are close for both linear and quadratic cases: o = 0.274 and ex = 0.267, 
respectively. From the diagram it is clear that errors for earthquakes with m > 6.5 deviate 
significantly from a linear trend. However, since the number of strong earthquakes is small, 
the residuals of the linear and quadratic cases do not differ significantly. 

The parameter values in ( ITT]) for two subsets of the CMT catalog are listed in Table [TJ 




(10) 




(11) 
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The e- values for earthquakes in the magnitude range 5.4-6.4 can be represented as 

e « 0.056 x i(T a54 ( m - 6 ) , (12) 

for shallow events in the 1982-2008 CMT catalog (see Tabled]). For small e the magnitude 
error a m is calculated as 

a m « . (13) 

1.5 log 10 V ; 

Modifying (jSJ) or Eq. 10 in Rhoades and Dowrick (2000), we obtain the following magnitude 
correction for the magnitude estimates perturbed by random errors: 

3 

m corr = m estim - -a^/3 1ogl0. (14) 

To apply this formula to b or correction, we estimate a m at two magnitude values (5.4 
and 6.4) and use ffl~4|) to compute 5 c m = m COTT — m estim . Performing such calculations for 
the relative moment error e in ( fi~3j) . we obtain the correction for (3 of shallow earthquakes 
0.0013: a (3 is decreased by about 0.2%. 

However, our calculations could not consider one important source of error. For many 
weak shallow earthquakes in the CMT catalog, no solution can be obtained for the tensor 
components M r g and M r ^ (Dziewonski et ai, 1981, p. 2829; Dziewonski and Woodhouse, 
1983; Frohlich and Davis, 1999). In such a case, E t q and E r § as well as M r g and M r< ^ are set 
to zero. About 4% of shallow earthquakes have this problem (Kagan, 2003, pp. 195-196). 
For strike-slip events which predominate in this group, the tensor components M r g and M r ^ 
are close to zero. This means that if the values of E^ and were available for these 
events in ffTU]) . the numerator value should be much greater, but the denominator would be 
essentially the same. This would significantly increase the resulting e-value. 

Moreover, apparently the relative moment error, e, is only part of the total seismic 
moment error. Dziewonski et al. (1981) and Dziewonski and Woodhouse (1983) suggested 
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the standard errors obtained in the CMT solutions likely to be underestimated. Kagan (2000, 
2002a) estimated that the reported errors are possibly 1/3 to 1/2 of the total. Because the 
bias in estimating (3 depends on the square of the magnitude estimation error (see Eq. [13]), 
a systematic bias as high as 1-3% may be caused by the decrease of the relative magnitude 
uncertainty with the increasing earthquake size. 

3.2 Earthquake sequences and their influence 

As mentioned in the Introduction, theoretical (3 estimates should be only relevant for earth- 
quake sequences, not individual events. We take that an earthquake belongs to sequences 
and sequences being the theoretical entity of interest, rather than individual earthquakes. 
Registration of aftershock sequences (Kagan, 2004; Enescu et ai, 2009) shows that immedi- 
ate aftershocks observed in high-frequency seismograms are included in mainshock or large 
aftershocks in catalogs based on long period registration. Thus, for example, the CMT 
catalog earthquakes include some close aftershocks. 

In our model earthquake sequences are produced by the same critical branching process. 
Later aftershocks are separated into individual events due to temporal delays controlled by 
Omori's law (Kagan and Knopoff, 1981). Occasionally, the first event in a sequence is weaker 
than the following events, in which case it is commonly called a 'foreshock' (ibid). Therefore, 
we could consider an earthquake cascade or foreshock-mainshock-aftershock sequence as one 
entity. 

Here we attempt to study the seismic moment distribution for earthquake sequences. To 
define the sequences, we use the results of the likelihood analysis of earthquake catalogs 
(Kagan, 1991b). We approximate an earthquake occurrence by a multidimensional Poisson 
branching cluster process. In this model spontaneous events are distributed according to 
the Poisson distribution, whereas dependent events in earthquake sequences are controlled 
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by a distribution characterized by a few adjustable parameters. Each event, even if it is 
a member of a cluster, may start its own sequence, hence aftershocks may be of the first-, 
second-, third-, etc., order. The parameters of the model are estimated through a maximum 
likelihood search. A similar scheme has been proposed recently by Zhuang et al. (2002, 
2004). 

As the result of the likelihood optimization, we evaluate probabilities (pij) of any i-th 
earthquake belonging to a j-th cluster or sequence (J2jPij — 1); Pa corresponds to the prob- 
ability that an earthquake will be considered independent. We use these probabilities (p^) 
to assign a part of the seismic moment of the i-th event to the j-th earthquake; the j-th 
earthquake might again belong to some k-th group, etc. This process continues until all 
earthquakes and their interconnections in a catalog have been counted. In the end, some 
of the aftershock moments are transferred to their mainshocks if the aftershock probably 
belongs to the particular mainshock's cluster. As a result of this seismic moment reassign- 
ment, some earthquakes may have a seismic moment below the magnitude threshold, m t . 
We remove these earthquakes from a catalog. 

As the number of earthquake sequences is always smaller than the number of earthquakes 
in a catalog, while the total moment in a catalog is constant, we should expect the f3- value 
for sequences to be smaller than that for individual earthquakes. For deep and intermediate 
earthquakes, the difference in the /3-values in calculations which use sequences and those 
using individual earthquakes is negligible. This small difference is due to a small number of 
aftershocks for these sequences in the CMT catalog (Kagan, 1999, Table 4). 

Generally, we can treat the probabilities of being independent (pa) as corresponding to the 
weight of an earthquake as it is included in calculations. However, to make our computations 
similar to those used for real catalogs, we simulate new catalogs leaving in only earthquakes 
whose pa exceeds a random number distributed uniformly in the [0-1] interval. Thus, we 
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obtain a 'declustered' catalog, in which we delete an earthquake according to its probability 
of being dependent. 

Table [2] shows several /3 measurements for two CMT shallow earthquake subcatalogs, 
1977-2008 and 1982-2008, with the magnitude threshold m t = 5.8 and m t = 5.6, respec- 
tively. For global datasets three types of computation were performed: (a) in the original 
list, (b) in a declustered catalog, where seismic moment has been preserved for each earth- 
quake, and (c) in a declustered catalog with aftershock moment transferred to an appropriate 
mainshock according to probabilities Pij (see above). We performed similar measurements 
for earthquakes in subduction zones (trenches) (Kagan et ai, 2010). Trench earthquakes 
have not been declustered, because some may have connections to outside events. Therefore 
the dependence probabilities can be biased. 

In Table [2] the 0- values are smaller by about 1-3% for the 1982-2008 dataset compared 
to the 1977-2008 catalog. Probably a higher average accuracy of these solutions (Kagan, 
2003) and larger magnitude range explains this reduction. Bird and Kagan (2004) showed 
that for global seismicity the minimum value of m c for some tectonic zones is of the order 
5.9-6.6. Extending the magnitude threshold to m t = 5.6 expands the power-law part of the 
plot, and the influence of the corner magnitude is smaller. 

The (3 bias in the Table is also caused by the mix of different earthquake populations 
with various corner magnitudes, m c (item 3 in the Introduction). This effect could explain 
why the /3-values for trenches are significantly lower (by about 5%) than the global ones. 
A global earthquake set consists of many populations, of which oceanic rift zones have the 
smallest m c - values (Bird and Kagan, 2004). These oceanic events are excluded from the 
subduction zone (trench) dataset. Hence, the estimate of (5 for trench earthquakes is closer 
to the theoretical value than for earthquakes in any other tectonic province. 

As expected, the (3- values decrease for the declustered catalogs, since excluded aftershocks 
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have smaller moment values. This reduction is even stronger for catalogs where the aftershock 
moment is assigned to their potential mainshocks. The /3-value decreases are about 4.5% 
and 8.5%, respectively. 

These bias estimates depend on the correctness of the calculations used to estimate 
earthquake probabilities. The likelihood procedure used to assign the probability of event 
independence is influenced by catalog quality, length and magnitude threshold. Given the 
presence of temporal boundaries, many relations between earthquakes are missed: some 
events of the beginning of the catalog may be aftershocks of previous strong quakes. Thus, 
instead of having a probability closer to zero, as they should, these aftershocks would have 
an independence probability equal to 1.0. 

Due to the magnitude threshold, some connections between events are not observable. 
Suppose there is a potential foreshock-mainshock pair: a larger earthquake is preceded by a 
smaller one. However, the first event is below the magnitude threshold and the larger quake 
is above it. Then this second event would be treated as independent; our calculations would 
not include this connection (Kagan, 1991b). 

Moreover, the likelihood model used in our inversion is not perfect (Kagan, 1991b). As a 
result, the independence probability values of earthquakes may not be fully counted leading 
to a bias in the j3 computations. Therefore, the reduction of the ,9-values due to the influence 
of aftershock sequences is likely to be greater. 

If we add up all the influences of aftershocks, we should see the /3-value decrease to about 
0.59-0.6. Kagan (1991a, p. 129) and Kagan (1999, Table 5) obtained a similar result for 
declustered shallow earthquakes. 
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3.3 Seismic moment tensor and its complexity 

The previous discussion assumed that scalar seismic moment is a fair measure of earthquake 
size. In reality seismograms are caused by excitation from many subevents during the main 
phase of an earthquake. Thus, the seismic moment tensor of an earthquake is a compound 
tensor sum of subevents. If all these subevents were identically oriented, the tensor sum 
would be proportional to a scalar sum of all the subevent scalar moments. However, detailed 
studies of earthquakes clearly indicate that subevent orientation changes significantly during 
ruptures. 

The Bulletin of the Seismological Society of America (BSSA) published several special 
issues dedicated to thorough analysis of several large earthquakes like the 1992 Landers, 
1999 Hector Mine, the 1999 Chi-Chi, Taiwan, the 2002 Denali, the 2004 Sumatra, and so 
on. These studies detail a very complex geometrical picture of the quake rupture process. 
The focal planes and slip vectors of earthquake subevents often rotate several degrees and 
even tens of degrees. Therefore, the seismic moment tensor solution and the resulting esti- 
mate of an earthquake scalar moment is subject to random fluctuations from the stochastic 
misalignment of earthquake components. 

In principle, we could avoid the systematic effect caused by source complexity if we used 
an earthquake's energy as a measure of its size. Energy is a positive scalar; thus, no bias due 
to source complexity would appear in the energy estimate. Unfortunately, estimates of the 
radiated seismic energy are highly uncertain and often differ by up to an order of magnitude 
(Perez-Campos et al, 2003). In contrast, the relative accuracy of evaluating seismic moment 
tensor is on the order of 10 a15 (Kagan, 2003). 

Scalar moment earthquake estimates should always be lower than the sum of the 
subevents' scalar moments. This occurs because of random fluctuations during earthquake 
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fracture. This effect would also bias upwards the estimated /3-values. Because we lack a 
comprehensive model of the earthquake rupture process which would enable us to estimate 
rigorously the resulting bias, we proceed by studying several approximations. These will 
give insight into the problem and provide an order of the magnitude estimate of possible 
systematic effects due to source complexity. 

We can estimate the influence of source complexity on the resulting estimate of earth- 
quake size by initially assuming that small elementary subevents have their sign selected 
randomly. If the sign changes with equal probability (p = 0.5), the resulting sum of the 
subevents would be an ordinary random walk. The walk would converge to Brownian mo- 
tion if the number of subevents, n, is large. The sum would be distributed according to the 
Gaussian distribution with the standard deviation proportional to \/n. The final value of 
the sum would be proportional to its standard error. (In a count of the Brownian sum, we 
use an absolute value of the final walk position. Therefore, the total 'moment' estimate is 
positive.) Therefore, in a critical branching process in which descendants are added with the 
random probability p = 0.5, the power-law index would increase by a factor of two: — 1.0. 

If we change the probability value from 0.5 to a higher level, p = 0.5 + 5, this would 
produce Brownian motion with a drift (Feller, 1966). For small n values, the walk behavior 
would resemble a regular Brownian motion, and later the sum would have a steady com- 
ponent n8. Thus, its behavior would be similar to the cumulative number increase with 
p = 1. 

Fig. [6] shows three simulated source-time functions with the time delay controlled by the 
Omori-type function. The cumulative functions for each curve are the sums of elementary 
subevents of the unit size and only the signs are different. The first with p = 1 (deterministic 
addition of events) is similar to Fig. 3 by Kagan and Knopoff (1981). The initial step-like 
increase of this function would likely be interpreted as a mainshock, whereas a few steps at 
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a later time would be labeled as aftershocks. The p = 0.75 function increases the same way 
as the first curve but with smaller amplitude; random fluctuations are not easily observable. 
The random walk function (p = 0.5) behaves more erratically and its total final ordinate is 
much smaller than that of the other two curves. Only the values of the curves at the end 
of a branching simulation, corresponding to the total moment of a sequence, are counted in 
our calculations. For the p-values close to 0.5, the random branching walk could end up as 
a negative cumulative sum; as we explained above, we take the absolute value of the final 
ordinate of a simulation run. These values are assumed to correspond to the total seismic 
moment of an aftershock sequence. 

Fig. [7] illustrates the above considerations. We simulated a critical branching process and 
counted the sum of events at the end of each simulation run, such as the extreme right-hand 
points of three curves in Fig. EJ These numbers are shown in the diagram in the log-log 
format. While the event numbers are small (less than 10), the discretization effects are 
noticeable. For the largest sequences, random fluctuations are observable, because there 
are few of these sequences. In the mid-number range, the deterministic number addition 
(p = 1) distribution (red solid curve) has an index (3 = 0.5. As expected, the Brownian 
walk addition (blue dashed curve) has (3 = 1.0. As explained above, the curves for the 
motion with a drift first follow the Brownian curve; then they are parallel to the red curve. 
Thus, in the beginning their index is (3 = 1.0, and for larger numbers it changes to (3 = 0.5 
confirming our predictions. The randomness in the number addition significantly increases 
the power-law distribution index. 

In Fig. [8] we show a more complicated test. In a critical branching process simulation, 
we sum up seismic moment tensors instead of scalar quantities. In each simulation run, we 
determine a norm of the total tensor sum which for the seismic moment tensor is equivalent 
to the scalar moment. Again the red solid curve shows the distribution when the tensors 
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are identical. For other curves the tensors are independently randomly rotated through the 
3-D rotation angle, $ (Kagan, 2003, 2009). The maximum angle for a double-couple focal 
mechanism is 120°; therefore, these tensors are rotated in a uniformly random manner. The 
angle $ = 0° corresponds to the zero rotation (red curve). If $ < 120°, then the rotation is 
restricted, being uniformly random only for angles smaller than $. 

The Fig. [8] diagram appears similar to the previous plot (Fig. [7]). If the tensors' orientation 
is identical, power-law exponent /3 = 0.5. For a completely random orientation, = 1.0, 
and for a restricted misalignment, the curves follow the latter distribution first and are then 
parallel to the former line. 

As discussed earlier, many earthquake ruptures exhibit significant variations in focal 
mechanisms. However, detailed analyses of individual earthquakes are still rare and insuffi- 
cient for rigorous statistical study. Therefore, we study the degree of misalignment in several 
mainshock/aftershock sequences. Kagan (2000) investigated the correlation of earthquake 
focal mechanisms and showed that the degree of mechanism 3-D rotation increases between 
earthquakes with temporal and spatial differences. Hence, we hope that immediate after- 
shocks of strong earthquakes will characterize the geometric complexity of their rupture 
process. 

To this end we studied all shallow (depth 0-70 km) earthquakes in the 1977-2008 CMT 
catalog with a magnitude ni\ = 7.5 and higher. All earthquakes (m > 5.6) are considered 
aftershocks within the first 7 days of m7.5 earthquake occurrence and closer than 

r = 75 x io( m i- 7 - 5 )/ 2 km, (15) 

(Kagan, 2002b) There are 105 mainshocks in the catalog and 81 of them have one or more 
aftershocks. 

To investigate the orientation differences between a mainshock and its aftershocks, we 
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calculate the correlation invariant or tensor dot-product J3 (Kagan, 1992; 2009) 

J 3 = E m ij n ij> ( 16 ) 

for the main event (ra^) and the sum of normalized tensors for the whole 7-days aftershock 
sequence (ny). Summation of repeating indices is assumed. Both and are normalized. 
In f|T6|) J 3 = 2.0 means that focal mechanisms are identical; J3 = —2.0 corresponds to 
components of both tensors having the opposite sign. 

A J3 histogram in Fig. [9] displays the correlation between tensors. Most correlation 
invariant values are close to 1.5-2.0. Thus, the aftershock focal mechanisms are similar 
to that of their mainshock. However, some ^-values are close to zero, and one is negative, 
testifying to a significant variation in the rupture process. The smallest J 3 -value is due to the 
November 2000 New Ireland earthquake sequence. The sequence started with a mw = 8.0 
left-lateral main shock on 16 November and was followed by a series of aftershocks with 
thrust mechanisms primarily (Geist and Parsons, 2005; Park and Mori, 2007). The negative 
J3- value signifies that the aftershocks have on average a slightly opposite orientation to their 
mainshock. 

Fig. [TU] displays two distributions of the ratio for the tensor sum of the mainshock and 
its aftershocks to the sum of their scalar moments 

R = I E m » I / E M , ( 17 ) 

where \rn,ij\ means the norm of the tensor. The aftershocks are selected according to the same 
criteria as in Fig. [9J If aftershocks have the same focal mechanism as the mainshock, the 
ratio would be 1.0. In the left diagram of Fig. [10] the moment tensors are not normalized; 
in the right plot they are. In the inversions of the earthquake rupture process (see the 
BSSA special issues mentioned above), several subevents of approximately equal size but 

23 



significantly different orientation in focal mechanism are often observed. This is the reason 
we investigate the normalized sums. 

The diagrams show significantly varied focal mechanisms in aftershock sequences. Large 
fluctuations are seen in the normalized sums especially. This result suggests again a con- 
spicuous randomness occurs in the focal mechanism orientation of earthquake sequences. 
By implication this should also occur during an earthquake rupture process. Such ran- 
dom fluctuations may noticeably decrease the measured earthquake size and influence the f3 
measurement. 

What is the size of the (3 measurement bias? All estimates shown above are indirect. 
Earthquake do not consist of an identically oriented or purely random collection of elementary 
sources. Various observations suggest that a rupture occurs over quasi-planar fault patches, 
so there will be a strong correlation between neighboring fault segments. This correlation 
is sometimes broken by significant fault branching. Kagan (1982) proposed a geometrical 
model of such stochastic rupture. Unfortunately, the degree of geometrical branching (0 O ) 
in this model is not well known for different tectonic provinces. Therefore, we cannot easily 
simulate and study such branching sequences. 

What /3-value change can be proposed as randomness result in the fault rupture orien- 
tation? Pure randomness yields /3- value increase by a factor of two (see Figs. [7] and [S]). 
Unfortunately, we cannot yet quantitatively study the complex geometry of earthquake rup- 
ture. We need to extrapolate from the measured misalignments of close aftershocks. These 
measurements indicate that complexity, though far from being completely random, is never- 
theless quite significant. For example, in Fig. [9] the correlation invariant is J3 = 1.46 ± 0.55. 
These values can be compared to purely random arrangements of double-couple focal mech- 
anisms (Kagan, 1992), where J3 = ± 0.89 has been obtained. 

The frequency plot for randomly oriented double-couples obtained by simulation is shown 
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in Fig. [TTJ As expected, the histogram curve is symmetric around J 3 = and reaches the 
maximum at abs(J 3 ) = 1. It would be interesting to obtain an analytical solution for J 3 , as 
we did for the 3-D rotation angle (Kagan, 2009). This can be done in the future. 

If focal mechanisms are all parallel, J 3 = 2.0 (see more in Kagan, 2009). Therefore, we 
can make a rough guess that the degree of the /3-value increase would be on the order of 
10-15%, when extrapolated to a time close to zero: during the mainshock rupture. This 
guess is obtained by comparing the average J 3 shift in Fig. [9] ( J 3 — 2.0 = 0.542) with that for 
the completely random arrangement in Fig. [TT] (J 3 — 2.0 = 0.0), and by comparing standard 
deviations for both cases: 0.55 and 0.89, respectively. 

Similar conclusions could be inferred from the results of the 3-D rotation angle distri- 
bution. The average angle between the mainshock focal mechanism and mechanisms of 
immediate aftershocks is on the order of 10° (Kagan, 2000). For a completely random ro- 
tation, the maximum angle is 120° and the average angle is 75.2° (Kagan, 2003). Given 
the source complexity, as demonstrated by the rotation angles, the (3 bias should be around 
8-12%. 

Analogous conclusions could be drawn from Fig. [T0J If immediate aftershocks had the 
same moment tensor solutions as their mainshock, the tensor/scalar sum ratio R in (JT7J) 
should be 1.0, and (1 — R) equal zero. Both average and standard deviation for (1 — R) in 
the plots display significant non-zero values. We infer that the source is complex and the /3 
bias may be on the order of few percent. 

To summarize the results of this subsection, we hypothesize that as a consequence of the 
random geometrical misalignment of a fault rupture, the measured /3-value may be increased 
by at least a few percent (3-7%) from its true size. The estimate above is conservative. More 
work needs to be done to obtain a more reliable value which could lead to an even greater 
/3-estimate decrease. 
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3.4 Centroid depth influence 

The CMT catalog supplies earthquake coordinates for the seismic moment centroid 
(Dziewonski et ai, 1981; 1983); the centroid is in the center of the moment release volume. 
The centroid distance from the fault edge cannot be smaller than a half of the earthquake 
rupture zone width. Closer to the surface or to the fault boundary the corner moment would 
approach zero. An inspection of seismic maps suggests that hypocenters of larger events 
are on average deeper than those for small earthquakes. Thus, the moment-frequency law 
(see Eqs. [5] and E]) would change due to an increase of the maximum earthquake size with 
depth (Kagan, 2002a, p. 539). Therefore the corner moment for deeper earthquakes would 
increase. 

As we explained in the Introduction, the depth accuracy for shallow earthquakes is 
presently insufficient to investigate observationally the dependence of the corner moment 
on depth. Therefore, we study a possible influence of the finite fault size by calculating 
a new distribution of earthquake size for a few simple models of the earthquake rupture 
pattern. These theoretical guesses would help evaluate the depth effect up to the order of 
magnitude. 

We assume that earthquakes are distributed over an infinite planar fault surface extend- 
ing either vertically for 20 km (imitating conditions for strike-slip faults in California), or 
distributed over an inclined fault with a width of 200 km (as in some subduction zones, see, 
for instance, Bird and Kagan, 2004). The variable corner moment for such faults is 

M' c = CO 3 , (18) 

where ( is distance from a centroid to a fault edge and C is an appropriate coefficient. 

In one model we assume that earthquake centroids are distributed uniformly over the 
fault surface. Then, using the algebraic and numerical facilities of mathematica (Wolfram, 
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1999), we calculate the new survivor function 

(19) 

where F(M) is defined by (El), L is a half-width of a fault plane, and T(. , .) is an incomplete 
gamma function (Abramowitz and Stegun, 1972, Eq. 6.5.3). 

The other possibility is to assume that earthquake centroid density increases linearly 
with increasing depth up to the middle of the fault width. The density decreases to zero 
thereafter. Kagan (2007, Fig. 6) shows that such a feature is a common occurrence. Then 
the survivor function would be 

h<Anr„ ( Mt Y t ( M t — M \ ( M t ^ £(-§, W) 

^(M) oc jF(M)Cd( = [-) J ^[-^^ td( = - . 

(20) 

Fig. [12] displays two survivor functions corresponding to Eqs. |6] and [201 In a loglog plot 
the former function has a linear part for the moment M values that are significantly smaller 
than the corner moment M c . We take a slope /3 to be 0.5. The curve has an exponential 
taper for M close to M c . On the other hand, the latter function ([2"U]) is slightly convex even 
for small moment values. It is formed by a sum of distributions similar to (J6j) but with the 
corner moment increasing from zero to the maximum, M c . Therefore, we can only calculate 
an effective slope (/?') of the curve; for the moment range 10 17 — 10 19 Nm, the slope j3' in the 
plot is 0.523. 

We calculate three theoretical curves for both equations f[T9l and 120]) : (a) fault width 
L = 200 km and M c = 10 23 Nm (m c = 9.33); (b) L = 200 km and M c = 10 22 Nm 
(m c = 8.67); and (c) L = 20 km and M c = 10 21 Nm (m c = 8.0). For the formula (EES]) 
the /^'-values are 1.3, 2.8, and 6.6% higher than the original /3-value, for the second formula 
these exponent increases are reduced by a factor of 1.4. 
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4 Discussion 



In the previous sections we analyzed the index of the power-law distribution for earthquake 
size (the (3- value) to argue that its true value is 1/2, the value suggested by theoretical 
arguments. The direct j3 measurements for scalar seismic moment based on catalog analyses 
(Kagan, 1999; Bird and Kagan, 2004) usually yield a value in the range 0.63-0.67, equivalent 
to the commonly known G-R 6-value of 0.95-1.0. Four systematic and random factors that 
bias the /3-value estimate upwards are investigated: dependence of errors on the magnitude, 
earthquake sequences, complexity of earthquake source, and a finite size of earthquake faults 
(items 4-7 in the Introduction). We found that these factors would decrease the observational 
/3-estimate by about 1-3%, 5-10%, 3-7%, and 2-6%, respectively. Of these values the third is 
most uncertain, because it is based on extrapolating immediate aftershock focal mechanisms 
to the mainshock' rupture time. 

If we combine the above biases and apply them to the most accurately determined (3- 
value, i.e., j3 equals from 0.63 to 0.64 for subduction zones (Bird and Kagan, 2004), the 
corrected /9-values would be on the order of 0.52-0.56. It is quite feasible that the second 
and third correction term are underestimated. This would imply that j3 is close to 0.5 and 
possibly equals \ exactly. 

What theoretical conclusions could be drawn from this result? Solid state physicists 
explain new results on the scale-invariant distribution of dislocation avalanche size by sug- 
gesting a new theoretical approach to crystal plasticity (Zaiser, 2006). According to this 
interpretation, at a micro-scale the crystal deformation proceeds through intermittent bursts 
similar to earthquakes. Only at a larger, meso-scale does plastic deformation proceed as a 
smooth, homogeneous, quasi-laminar flow process. Crystal boundaries seem to influence this 
transition. In a single crystal the power-law distribution for energy of dislocation avalanches 
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is observed at the scale range of 10 6 . In polycrystal materials, the power-law distribution of 
bursts is also observed, but its size is limited by an upper cutoff. 

Zaiser and Moretti (2005) and Csikor et al. (2007) propose the following probability 
density function for the dislocation avalanche energy or strain 



where C is a normalization constant, r is a scaling exponent (r ~ 1.5), and sq is the char- 
acteristic strain of the largest avalanches. This formula is similar to our (JHJ) with exponent 
t = 1 + (3, but the decay taper at large strain values is Gaussian like. We use instead the 
exponential decay. Because statistics on the largest events are insufficient in both cases, we 
cannot distinguish by observation between these formulas. 

Therefore, plastic, ductile deformation proceeds by two very different mechanisms: (a) 
intermittent displacement at micro-scale with scale-invariant distribution of strain steps and 
the universal value of the power-law exponent (r = 1.5), and (b) in contrast a smooth flow 
at larger scales. Because detailed quantitative observation at small sub-grain scales was not 
possible until recently, the first mechanism had been largely ignored (Zaiser, 2006, p. 241). 

The above considerations can be supported to some degree by recent analysis of earth- 
quake size distribution. Bird and Kagan (2004) found that the exponent of the power-law 
distribution appears universal in all eight tectonic provinces of global seismicity. In stark 
contrast, the corner moment differs by many orders of magnitude, from 10 18 (m c = 6.0) for 
oceanic normal faults to 10 23 ' 3 Nm (m c = 9.5) for subduction zones. 

Kagan (2002a, pp. 538-9) proposed that the observed b- value differences in volcanic areas, 
at creeping faults, and at the boundary between brittle crust and plastic deformation in the 
upper mantle may also be due to significantly varied corner moments. If earthquake popula- 
tions with different m c are mixed, the resulting statistical distribution could be interpreted 




(21) 
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as belonging to a power-law with the exponent /3 significantly exceeding 0.6-0.7, the value 
normally observed in tectonic earthquakes. 

Another geophysical phenomenon, the non-volcanic seismic tremor (Schwartz and 
Rokosky, 2007; Beroza and Ide, 2009, and references therein), may be explained by the 
same physical mechanism. The tremor represents long duration (minutes to hours) of a 
high-amplitude seismic signal, appearing similar to many small concatenated earthquake 
signals (Shelly et al., 2007). The first observation of such non- volcanic tremors came during 
the aftershock sequence of the April 26, 1966 Tashkent (Uzbekistan) earthquake (Antsyferov 
et al, 1971a, 1971b). 

The tremor signals are sometimes quite pulsed in nature. For example, the temporal 
cumulative plot of seismic moment increase for tremors (Fig. 4d in Hiramatsu et al, 2008) 
look similar to curves of crystal micro-deformation due to dislocation avalanches (Fig. 14 in 
Zaiser, 2006): both diagrams show that the displacement increases in discrete steps, each 
step is followed by a plateau. 

Such tremors have been registered in diverse tectonic environments recently (Japan, 
Cascadia, New Zealand, Costa Rica, Taiwan, California). Tremor and other slow-slip events 
are typically found on the deep extension of faults, just below the region that produces the 
more familiar, 'ordinary' earthquakes. This recent observation of tremors has resulted in a 
flurry of research across many geologic and geophysical disciplines. 

If tremors are a feature transitional between real earthquakes and seismic signal bursts, 
they are similar to the dislocation avalanches described above. Both phenomena occur in 
conversion from the brittle to plastic mode of solid deformation. Seismic tremors, which 
are interpreted as small, continually occurring earthquakes, may also have the same scale- 
invariant, power-law features as earthquakes in brittle crust. Hiramatsu et al. (2008) mea- 
sured moment-frequency relation for tremors and found that it can be approximated by an 
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exponential distribution rather than a power-law. However, because of low signal-to-noise 
ratio for tremors, only the upper tail of tremor size distribution can be observed, and the 
upper tail for earthquake size distribution (see Eq. [5]) is also exponential (Kagan, 2002a; 
Bird and Kagan, 2004). The higher end of the size distribution for dislocation avalanches 
also exhibits a non-power-law dependence ( J2T1) which is close to exponential law. We may 
conjecture that as with dislocation avalanches, the size distribution of smaller tremor events 
would be the power-law with the universal value of the exponent constant ((3 = 0.5). Further 
study of tremors should answer this question. 

If we are correct about universality of the 0- value constant (/3 = 0.5), the observed vari- 
ations in the 6-parameter result from systematic and random effects not properly accounted 
for (see the Introduction, items 1-3). Therefore, all attempts to connect 0- value variabil- 
ity with various physical parameters and conditions are eventually bound to fail. However, 
studying the b- or /3-values in local and regional earthquake catalogs may still be useful, es- 
pecially if such investigations are needed to evaluate seismic hazard and seismicity forecasts 
that would be prospectively tested with the same catalogs. In addition, when seismic activity 
or earthquake productivity level is calculated for large earthquakes, the regular /3-estimates 
can be used. 

If the hypothesis about the (3- value constancy is correct, we should investigate spatial 
changes in the corner moment that seem to explain major modifications in the deformation 
processes in solids. Bird et al. (2009) showed that the relation between relative plate velocity 
and seismicity is non-linear for several types of plate boundaries. Can the change in corner 
moment explain some of these non-linearities? 

If the hypothesis that the power-law exponent is a universal constant and the corner 
moment is variable is correct, then it would provide a new theoretical approach to features 
of earthquake occurrence and account for the transition from brittle to plastic deformation. 
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More extensive investigation of corner moment behavior may afford new insight into regular 
earthquake occurrence and recently discovered slow deformation and seismic tremor events 
at the brittle-plastic crust boundary. As often happens in complex systems, new laws and 
features may be found to illuminate the transition from brittle fracture to plastic flow. 
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Table 1: Parameter values for relative errors. 



Approximation n Co c\ c 2 pa 



CMT 1977-2008, m t = 5.6 
Linear 8508 -1.20 -0.352 -0.47 0.274 

Quadratic 8508 -1.23 -0.513 0.198 -0.47 0.267 



CMT 1982-2008, m t = 5.4 
Linear 11600 -1.21 -0.447 -0.57 0.274 

Quadratic 11600 -1.26 -0.537 0.208 -0.57 0.266 



n, the number of m > m t events; for c , Ci, and c 2 see ffTTl) ; 

p, correlation coefficient; a, standard deviation of fit. 

For the linear approximation, the c 2 parameter is shown as a dash. 
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Table 2: The /3-values for various subdivisions of CMT catalog. 



# 


Earthquakes 


m t = 


5.6 


m t = 


5.8 






P 


Eq. # 


P 


Eq. # 


1 


2 


3 


4 


5 


6 


I. 


Global 


0.6773 


7369 


0.6820 


5450 


2. 


Global declustered 


0.6480 


5841 


0.6568 


4498 


3. 


Global aftershocks included 


0.6229 


5605 


0.6366 


4358 


4. 


Trenches (Subduction zones) 


0.6463 


4805 


0.6507 


3223 



Shallow earthquakes in 1977-2008 (m t = 5.8) and 1982-2008 (m t = 5.6) CMT catalog. 
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Figure 1: 

Number of earthquakes with moment (M) larger than or equal to M as a function of M 
for the shallow earthquakes in the CMT catalog during 1977-2008, moment threshold M t = 
10 17 ' 7 Nm (m t = 5.8), the total number of events 5450. Power-law approximation (equivalent 
to Gutenberg-Richter law) is shown by dotted line. Dashed line shows tapered Gutenberg- 
Richter distribution: the G-R law restricted at large seismic moments by an exponential taper 
with the corner magnitude m c = 8.9. The slope of the linear part of the curve corresponds 
to = 0.68. 
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Figure 2: 



Log-likelihood maps for the distribution of scalar seismic moment of earthquakes: The CMT 
catalog time span is 1977 January 1 - 2008 December 31; the seismic moment cutoff is 10 17 7 
Nm {m t = 5.8); the number of events is 5450. iJ-sign on the plot denotes the maximum 
likelihood estimate of the parameters of interest. 
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Figure 3: 

Difference between observed magnitude-frequency relation and its approximation by tapered 
G-R law (see Fig. [1]). CMT catalog 1977-2008, magnitude threshold m t = 5.8. The dashed 
line indicates where the difference is based on fewer than 10 events. 
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Fig. 4 
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Figure 4: 



Percentage of moment tensor solutions inversion of CMT data based on different waves: 
body waves - red line, circles; surface waves - blue line, x-marks; mantle waves - green 
line, pluses. Dashed lines are for data with fewer than 10 earthquakes in the 0.1 magnitude 
interval. 
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Figure 5: 



Relative error e versus moment magnitude for shallow earthquakes m t = 5.6 in the 1977-200^ 
CMT catalog. The curves show two approximations: linear and quadratic fits. 
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Fig. 6 
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Figure 6: 

Simulated source-time functions in a critical branching process. Red dotted line - positive 
number addition (p = 1.0). Green dashed line - unequal positive/negative number addition 
(random walk with a drift, p = 0.75). Blue solid line - equal positive /negative number 
addition p = 0.5 (random walk). 
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Fig. 7 
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Figure 7: 



Distribution of event numbers in a critical branching process. Red line, circles - positive de- 
terministic number addition (p = 1.0). Blue line, x-marks - equal positive/negative number 
addition p = 0.5 (random walk). Green line, pluses - unequal positive/negative number ad- 
dition (random walk with a drift, p = 0.51). Cyan lines, squares - unequal positive/negative 
number addition (random walk with a drift, p = 0.55). 
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Figure 8: 

Frequency plot of tensor sum norm in a critical branching process. Red line, circles - tensor 
sum with no rotation, the 3-D rotation angle, $ = 0°. Blue dashed line - random rotation 
($ = 120°). Green line, pluses - limited random rotation, $ = 80°. Cyan line, squares - 
limited random rotation, $ = 30°. 



49 



Fig. 9 
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Figure 9: 



1.5 



Frequency plot of correlation tensor invariant (tensor dot-product) for m7.5 mainshocks and 
sum of their immediate aftershocks in 1977-2008 CMT catalog. Average J 3 = 1.458, its 
standard error oj = 0.5535. 
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Fig. 10 
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Figure 10: 

Frequency plot of tensor/scalar sum ratio for m7.5 mainshocks and immediate aftershocks 



in 1977-2008 CMT catalog: (a) Unnormalized sum - average (1 — R) — 0.0128, its standard 



error cr = 0.0627. (b) Normalized sum - average (1 — R) = 0.1642, its standard error 
a = 0.1639. 
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Figure 11: 

Frequency plot of tensor dot-product invariant for random rotation of double-couple sources. 
Average J 3 = 0, its standard error aj = 0.8945. 
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Figure 12: 

Two theoretical moment-frequency curves: the number of earthquakes with moment (M) 
larger than or equal to M as a function of M, moment threshold M t = 10 17 Nm (m t = 5.33). 
Dashed line shows tapered Gutenberg-Richter distribution: the G-R law restricted at large 
seismic moments by an exponential taper with the corner moment 10 210 Nm (m c = 8.0). 
The slope of the linear part of the curve corresponds to /3 — 0.50. The solid line is a plot for 
Eq. [19j with a half-width of a fault, L = 10 km, and C = 1.0. The curves are normalized so 
that the solid line has an ordinate 1.0 at the right-hand end. 
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